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l> . ABSTRACT 

O 
O 

The fragmentation process of primordial-gas cores during prestellar collapse 
is studied using three-dimensional nested-grid hydrodynamics. Starting from the 
^ | initial central number density of n c ~ 10 3 cm -3 , we follow the evolution of rotat- 

ing spherical cores up to the stellar density n c ~ 10 22 cm -3 . An initial condition of 
the cores is specified by three parameters: the ratios of the rotation and thermal 
energies to the gravitational energy (/3 , and a , respectively), and the ampli- 
tude of the bar-mode density perturbation (A^). Cores with rotation /3 > 1CT 6 



Oh 



O 

are found to fragment during the collapse. The fragmentation condition hardly 

depends on either the initial thermal energy «o or amplitude of bar-mode pertur- 
bation A^. Since the critical rotation parameter for fragmentation is lower than 
that expected in the first star formation, binaries or multiples are also common 
Q\ ' for the first stars. 

O 

Subject headings: cosmology: theory — galaxies: formation — hydrodynamics- 
stars: formation 

O 
> 

More stars are in binaries or multiples than in singles in the solar vicinity (Duquennoy & 
Mayor 1992; Fischer & Marcy 1992). Is this also true of the first-generation stars? Being very 
massive (> 1OOM ; see e.g., Bromm & Larson 2004), first-star binaries, if they are formed, 
can contribute greatly as progenitors of gamma-ray bursts (Bromm & Loeb 2006; but see 
Belczynski et al. 2006 for the opposite view) or gravitational wave sources (Belczynski 
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et al.2004). However, binary formation in primordial gas has sometimes been considered 
unlikely, based on the fact that no cooling mechanism induces fragmentation of high-density 
(> 10 4 cm~ 3 ) cores (e.g., Ripamonti & Abel 2004). 

Rotation can change this picture. Unlike the case of present-day star formation, where 



the a ngular mom entum is effecti vely transported by magnetic braking (IBasu fc Mouschovias 



19941 ) or outflow (ITomisa ka 200 2[), in first-star f ormation, the magnetic field is extremely weak 



(e.g., I'l'akahashi et al.M2005l ; llchiki et al.l 120061 ) and thus expected to play little role in the 
angular momentum transport. With conserved angular momentum, the effect of rotation is 
enhanced in the course of the collapse of the core, and a disk is expected to form. Do such 
disks fragment to binaries or multiples? To answer this question, by three-dimensional hy- 
drodynamics, Saigo et al.(2004) demonstrated that primordial-gas cores with rapid rotation 
indeed fragment into binaries. Nonetheless, their analysis is not satisfactory in the follow- 
ing respects: First, they started calculation at rather high density ~ 10 10 cm -3 , and thus 
were forced to assume somewhat ad hoc initial conditions of high angular-momentum cores. 
In practice, the angular momentum might have been transferred by the non-axisymmetric 
structure before reaching this density. In addition, they presented only one model of a core 
experiencing fragmentation. Recently, Clark, Glover, & Klessen (2007) performed three- 
dimensional hydrodynamical simulation of star formation in a low-metallicity turbulent 
medium and found that rotation causes fragmentation even in metal-free gas, as suggested 
by Saigo et al. (2004). However, they do not discuss fragmentation conditions of the cores 
in their analysis. 

In this paper, we calculated the evolution of primordial-gas rotating cores for twenty 
orders of magnitude in the density contrast, starting from the formation of dense cores 
(n ~ 10 3 cm -3 ) up to stellar-core formation (n ~ 10 22 cm -3 ). We adopt a realistic equation 
of state derived by Omukai et al. (2005), rather than a simple polytropic relation P oc p 11 
used by Saigo et al. (2004). To determine the fragmentation conditions for the primordial- 
gas cores, we calculated 38 models with different values of rotation, thermal energy, and 
amplitude of non-axisymmetric perturbation. Using the initial ratio of rotation to gravi- 
tational energies /3q, the conditions for a primordial-gas core to fragment are /3q > 10~ 6 , 
while dependence of these conditions on either the thermal energy or the amplitude of non- 
axisymmetric perturbation is very weak. These conditions are easily satisfied in the first-star 
forming cores, and thus binary formation is also important for first stars. 

The following is the plan of this paper. In this paper, after first summarizing the basics 
of our model in §2, we describe the numerical method in more detail in §3. We present the 
results for the initially quasi-stable cores in §4, and those for highly gravitationally unstable 
cores in §5. The fragmentation criterion is examined in §6. Finally in §7, summary and 
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discussions are presented. 

2. Model Settings 

We solve the equations of hydrodynamics including self-gravity: 

f t +V-(pv) = 0, (1) 

P -^ + p(vV)v = -VP-pV<j ) , (2) 

V 2 = 4vrGp, (3) 

where p, v, P, and denote the density, velocity, pressure, and gravitational potential, 
respectively. For gas pressure, we use a barotropic relation that approximates the result 
for a primordial-gas core by Omukai et al. (2005). Omukai et al. (2005) studied the 
thermal evolution of low-metallicity star-forming cores using a simple one-zone model for 
the dynamics, where the core collapses approximately at the free-fall rate and the size of the 
core is about the Jeans length, as in the Larson- Penston self-similar solution (Penston 1969; 
Larson 1969). On the other hand, detailed thermal and chemical processes in primordial gas 
are treated. Their result and our fit to it are plotted against the number density in Figure [TJ 
To emphasize variations of pressure with density, we plot P/n, which is proportional to the 
gas temperature divided by the mean molecular weight. 

As an initial condition of the cores, we use a density profile increased by a factor of 
/(> 1) from that of the critical Bonnor-Ebert (BE) sphere (Ebert 1955; Bonnor 1956). For 
the critical BE sphere, the central number density is n Cj o = 1 x 10 3 cm -3 , and the initial 
temperature is about 250K at the center from our adopted barotropic relation (Fig{T]). The 
critical radius of the BE sphere is R c = 6.45 c s /^AtiGpbe(0) = 6.5 pc. Outside this radius, 
uniform gas density of 71be(-Rc) = 71.1/ cm -3 is assumed. The total mass contained in the 
critical BE sphere is M c = 6.2 x 1O 3 M , and thus our initial core is / times more massive 
than this value. In Figure [2J the profiles of density and cumulative mass of the initial state 
are plotted. 

Initially the core rotates rigidly with angular velocity Qq around the z-axis. The initial 
models are characterized by three non-dimensional parameters: the ratio of the thermal to 
gravitational energy («o), the ratio of the rotation to gravitational energy ((3q), and the 
amplitude of the non-axisymmetric perturbation (A^). The models with different (a , (3q) 
are constructed by changing the initial density enhancement factor (/) and angular velocity 
(Qq). Parameters «o an d A) are summarized in Tabled], along with the density enhancement 
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factor /. Non-axisymmetric density perturbation of the m = 2 mode, i.e., bar mode, is 
added to the initial cores. In summary, the density profile of the core is denoted as 



p(r) 



PbeM (1 + Sp) f for r < R c 
p BE (R c )(l + 5p)f for r>R c 



(4) 



where Pbe(^) is the density distribution of the critical BE sphere, Sp is the non-axisymmetric 
density perturbation. For the bar(m = 2)-mode, 



5p = A^ir I ' R c ) cos 20 for r > R c 



(5) 



where A^ is the amplitude of the perturbation. We use A^ as a free parameter, and study 
cases with A$ = 0.001 — 0.3 (see Table HI). To distur b the m = 2 symmetry in a high-density 
region (for detail, see iMatsumoto &: Hanawal 120031 ) . we add the m = 3 mode of velocity 
perturbation as 

tt = tt [1 + tt 3 (r/R c ) cos 30] . (6) 

The amplitude ^3 is set to 10~ 5 in all the models. 



Numerical Method 



We adopt a nested grid method (for details, see Machida et al. 2005a; 2006a) to obtain 
high spatial resolution near the center. Each level of a rectangular grid has the same number 
of cells (= 256 x 256 x 16), with the cell width h(l) which is a function of grid level I. The 
cell width is halved at each increment of the grid level. The highest level of grids increases 
with time: a new finer grid is generated whenever the minimum local Jeans length Aj falls 
below 8 h(l majK ), where /i(/ max ) is the cell width of the current finest grid. The maximum level 
of grids is restricted to Z max < 30. The generation of a new grid at the density maximum 
ensures the Jeans condition of iTruelove et al.l (119971 ) with a margin-of-safety factor of two. 
We begin our calculations with four grid levels (i.e., Z max = 4). The box size of the initial 
finest grid I = 4 is chosen to be 2R C , where R c is the radius of the critical BE sphere. The 
coarsest grid (I = 1) thus has a box size of 2 4 i? c . The mirror symmetry with respect to the 
z=0 is imposed. At the bounda ry r = 2 4 R r , the ambient gas is set to rotate with angular 
velocity Qq (for more detail, see IMatsumoto fc Tomisaka 1 12004 ) . 



4. Evolution of quasi-hydrostatic cores 

According to numerical simulations (Bromm et al. 1999; Abel et al. 2002; Yoshida et al. 
2006), in the case of first star formation, self-gravitating isolated cores form at the number 



-5- 



density 10 3 — 10 4 cm -3 , where temperature takes its minimum (see Fig. [I]). These cores are 
close to the hydrostatic equilibrium, i.e., a — 1> which is a general property of cores forming 
by gravitational instability. In this section, we study such cases as ao — 1, and we defer the 
effect of different a to Sec. 5. 

In Figure [3j the final states of our simulation for the cores with «o — 1 are presented for 
various combinations of initial rotation parameter f3 and bar-mode perturbation A^. The 
model parameters (a , Po, A$) are summarized in Tabled] (models 1-24 for those in Fig. [3]). 
Also presented in Table [TJ are the number density (rif) at fragmentation, the number of 
fragments, the separation between the fragments at the end of the calculation, and the mode 
of fragmentation (i.e., ring or bar-mode; see below). Among models presented in Figure [31 
for models 1-4 with the highest rotation parameter (/5 = 0.1), we need to adopt values of 
a (= 0.7) smaller than in the other models (a = 0.83) to initiate the collapse. In Figure [3] 
we see that cores with rotation (3q > 10~ 5 fragment for all values of bar-mode perturbation 
Acf, in our calculated range. 

It is known that cores fragment only after they become either thin disks or elongated 
filaments (Tohline 1980). Correspondingly, the mode of fragmentation can be classified into 
two types (See Fig. H|). One is ring- mode fragmentation induced by rotation, and the other is 
bar- mode fragmentation induced by bar instability (Machida et al. 2004; 2005b). The mode 
of fragmentation is determined by the axis ratio. If the axis ratio does not grow sufficiently 
during the collapse, the core becomes ring-like due to the rotation and then fragments into 
pieces (i.e., ring-mode fragmentation) . On the other hand, a core with small rotation but with 
a large enough bar-mode perturbation results in an elongated structure, and fragments due to 
the bar instability (i.e., bar-mode fragmentation) . For example, models in the second column 
of Figure [3] experience ring- mode fragmentation, while in the fourth column where /3q > 10~ 3 
the cores fragment through the bar mode. The mode of fragmentation determines the fate 
of the fragments to some extent. Survival of fragments and formation of binary/multiple 
are easier in the ring mode due to large orbital angular momentum between the fragments, 
while in the bar mode, they tend to merge, and fragme ntation is repeated in the merged 



core (IMatsumoto fc Hanawal 120031 ; iMachida et al.ll2005bl ). 



Our calculations are terminated when either the Jeans condition is violated in the max- 
imum (/ max = 30) grid level or the maximum density exceeds n c > 10 22 cm~ 3 . Since the 
former situation easily occurs when the fragments leave the center, where the grid is finest, 
after the fragmentation we can continue calculation only for a short duration. Due to this 
limitation, our calculation covers only a rather low-density range for models with high /3 
(models 1-4), where fragmentation occurs at low densities (n c ~ 10 11 cm~ 3 ). In these mod- 
els, the cores might experience another episode of fragmentation if we continued calculation 
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further. For the other models in Figure [31 our calculation covers the entire evolution until 
the stellar core formation after the equation of state becomes adiabatic. 



4.1. Case Studies 

In this section, we see the evolution of cores in models 22, 6 and 8, which are exam- 
ples of a non-fragmenting core, a core fragmenting in the ring-mode, and in the bar-mode, 
respectively. 

We define here the following three quantities which describe the state of cores at a 
given time for later use. To examine the evolution of rotation velocity around the center, we 
use angular velocity normalized by the free-fall timescale; u c = fi c /(47rGp c ) 1 / 2 . Hereafter, 
we call this quantity the normalized angular velocity. For a rigidly rotating sphere with 
uniform density, the normalized angular velocity is related to the ratio of the rotation to the 
gravitational energy (j3 c ) by 

tf c R 3 nl 

^ c 3GM 4nG Pc ' [) 

where R and M are the radius and mass of the sphere, respectively. Even when the core 
deviates somewhat from the sphere in the course of the collapse, equation (|7|) remains a good 
indicator for rotation (Machida et al. 2005a; 2006a). Also, to see the degree of deformation 
from a sphere and approach to a disk, we define oblateness of the core as e b = {hih s ) l l 2 /h z , 
where hi, h s , and h z are, respectively, lengths of the major, minor, and z-axes derived from 
the moment of inertia for the high-density gas p > 0.1p c , i.e., density higher than 10% of 
the central value. Also, as the degree of deviation from a circle when projected on the x — y 
plane, the axis ratio is defined as e ar = hi/h s — 1. 



4-1.1. Non-fragmenting Core: Model 22 

As an example of cores that do not experience fragmentation, we describe the evolution 
in model 22, whose parameters are (ao, Po, A^) = (0.83, 10~ 6 , 10 -2 ). In this case, a 
hydrostatic core ("stellar core") forms without fragmentation. This structure is shown in 
several evolutionary stages in Figure Panel a shows the initial state. We adopt the critical 
BE sphere as an initial state (Sec. 2), but add 1% of density perturbation for triggering 
gravitational collapse. The initial density contrast is thus 14 between the center and the 
ambient medium. 

Since the rotation is very slow (/3q = 10~ 6 ) in this model, the collapse becomes almost 
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spherically symmetric in the early phase, as seen in panels b-d. In fact, the infall motion is 
almost radial, i.e., the radial component of the gas velocity is much larger than the azimuthal 
component, v T 3> v^. The effect of rotation becomes important with the collapse, and 
deviation from the spherical collapse becomes obvious around n c ~ 10 17 cm -3 (see Figs. E]e, 
where n c = 7.1 x 10 18 cm -3 ). Now the azimuthal velocity is comparable to the radial one 
(t>0 ~ v t ), and the gas falls to the center with significant rotation. 

In Figure El we plot the normalized angular velocity (Fig. Eh), oblateness (Fig. |B]&), 
and axis ratio (Fig. [B]c) against the central number density for this case (model 22) along 
with those with different /3q but the same A^ = 10 -2 (models 2, 6, 18) for comparison. 
Those quantities are calculated for the central region whose density is higher than 10% of 
the central value. 

For the spherical rigid-body rotation, angular velocity increases as Q, c oc nj owing to 
the conservation of angular momentum, and thus normalized angular velocity to c oc n c ■ In 
model 22, this spherical collapse phase continues until n c ~ 10 19 cm -3 , where the saturation 
of the normalized angular momentum around u c ~ 0.2 — 0.3 is observed. This indicates that 
a thin disk begins to form at the center, for which the normalized angular velocity increases 
as uj c oc nZ , which is almost constant in our case where 7 ~ 1.1. This can be derived 
from the conservation of specific angular momentum j/M = irzu 2 Q c /7r'uu 2 'E = const., where 
w is the radius of the disk and the column density £ can be expressed as X ~ pH using the 
disk scale height H ~ c a ) 'y/Gp oc p 7-1 ^ 2 , where c s is the sound speed. 

The formation of the disk can be seen in the behavior of the oblateness. For n c < 
10 19 cm -3 , the collapse is spherical and the oblateness remains close to unity (see Fig. [6]6). 
After the normalized angular momentum reaches the saturation level and the radial collapse 
is hindered to some extent by centrifugal force, the oblateness begins to increase. When the 
central density reaches n c = 3xl0 21 cm -3 , a stellar core with an initial mass of 6.7 x 1O~ 3 M , 
forms at the center, surrounded by a shock. Figure [5jf shows the structure immediately after 
the stellar core formation. In this figure, the surface of the stellar core is indicated by the 
thick black line. The density and mass of the stellar core at the formation in our model are 
well in agreement with the results of Omukai & Nishi (1998), who calculated the collapse of 
a primordial star-forming core including detailed radiative and chemical processes under the 
assumption of spherical symmetry. After the formation, the stellar core becomes round in 
shape and the oblateness stops increasing. The oblateness can grow remarkably only after 
the saturation of uo c and before the formation of the stellar core. In this model, the saturation 
of uj c occurs only slightly before the stellar core formation, owing to a small initial rotation 
(3q. Thus, the core does not become oblate enough to fragment in the ring mode (see Sec 
4.1.2). 



Next, we see the growth of the bar-mode perturbation e ar in Figure [6]c. Recall that we 
add 1% of bar-mode density perturbation to the initial state, i.e., A^ = 10~ 2 . The axis ratio 
of the core in model 22 remains very small throughout the evolution (e ar < 10~ 2 ). The growth 
of the bar-mode perturbat i on ca n be described by a linear perturbation theory developed 



by lHanawa fc Matsumoto I (120001 ) and lLail (120001 ) for small deviation from the spherical self- 



similar solution. They showed that a collapsing spherical core with a polytropic index 7 
is unstable against a bar-mode perturbation when 7 < 1.097, and the perturbation of the 
bar mode grows as e ar oc n s c , with s being a constant for a given 7. For the isothermal case 
(7 = 1), for example, s = 0.177 (^ 1/6): the growth rate is the same as that of the normalized 



angular velocity uj c (lHanawa &: Matsumoto Ill999l ). For more reference, see Figure 14 of ? 



where the power index s is shown as a function of 7. In our case, the polytropic index ranges 
in 1 < 7 < 1.15 (see Fig. [[]), corresponding to the range of the power index —0.4 < s < 0.2. 
For 10 8 cm -3 < n c < 10 16 cm -3 , 7 becomes ~ 1, and then e ar grows (Fig. [H]c), although the 
growth rate is somewhat lower than in the isothermal case. 

We show the evolution of the radial density distribution averaged around the rotation 
axis in Figure Ha. With only slow rotation, the evolution of both density and velocity 
distributions resembles closely a result of the spherically symmetric calculation (Omukai & 
Nishi 1998). The collapse proceeds in a self-similar way with a polytrope index 7 ~ 1.1: the 
evolution of density distribution resembles closely the Larson-Penston self-similar solution, 
whose density gradient in the envelope is oc r~ 2 ' ( 2 " ~ r~ 2,2 . The appearance of a shock front 
due to the stellar core formation is seen in the velocity distribution at state 6 (Figure [7]o). 



4-1.2. Ring-mode Fragmentation: Model 6 

Here, we see an example of a core with rapid rotation. In Figure [8], we show the 
evolutionary sequence of model 6, whose parameters are (ao, Po, A^) = (0.83, 10~ 2 , 10 -2 ). 
Rotation energy of this model is 10 4 times larger than that of model 22 above, while the 
other parameters A^ and a^ are the same. In contrast to model 22, where the collapse is 
almost spherical, here rotation plays a major role. After some contraction, the core becomes 
ring-like, owing to the centrifugal support, and eventually fragments. 

The core is still spherical early in the evolution (Fig. [8]a; e b = 1.2, n c = 3 x 10 7 cm~ 3 ). 
With further contraction, rotation makes it oblate perpendicular to the z-axis (Fig. [S]&; £ b — 
3.8, n c = 2x 10 11 cm -3 ). As in model 22 above, the normalized angular velocity grows initially 
as uo c oc n c and then saturates around oo c ~ 0.2 — 0.3 (Fig. [6]a). In this model, however, 
the saturation occurs earlier at n c ~ 10 7 cm~ 3 than in model 22 (~ 10 19 cm~ 3 ), owing to 
faster initial rotation. The normalized angular velocity oscillates around the saturation value 
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uj c ~ 0.2 — 0.3 thereafter, which corresponds to j3 c = 0.04 — 0.09 by the relation ([7j). In other 
words, after the ratio (3 C reaches a certain level of 0.04 — 0.09, the core contracts, keeping a 
constant ratio between the centrifugal and gravitational forces near the center. The approach 
to a disk-like configuration can be observed in the increase of oblateness (the dashed line in 
Fig. [6]6) after the saturation of the normalized angular velocity. The oblateness also saturates 
at e b — 10 and then oscillates after that. Note that the face-on view on the x — y plane 
remains quite circular even after the thin disk formation. For example, the axis ratio e aT is 
~ 10~ 3 at the time of Panels a and b of Figure El 

As can be seen in the middle and lower panels of Figure [8]c, two disks are nested, and 
each disk is sandwiched by horizontal shocks above and below. The outer and inner shocks 
exist at z — ±0.2 AU and z = ±0.05 AU, respectively. Outside the outer disk, cavities 
exist along the z-axis. The emergence of the shocks is related to the oscillatory behavior of 
oblateness after its saturation: a pair of shocks appear whenever the core becomes maxim ally 



oblate and then bounces back (IMatsumoto et al.l 119971 ; iMachida et al.l l2005al . l2006al ). In 



this model, two maxima of the oblateness appear around n c ~ 10 13 cm -3 and 10 17 cm -3 
(Fig. \5$)). Correspondingly, the outer shock appears at n c ~ 10 13 cm -3 and the inner shock 
at n r — 10 17 cm~ 3 . 



<c 



The axis ratio of an oblate core grows faster than the linear theory predicts (see Fig. [6]c). 
Although the core deforms to an ellipsoid by this effect (Fig. [8]c; n c ~ 5 x 10 18 cm^ 3 ), 
the deformation is not sufficient owing to the initial small value of A^. Consequently, the 
core eventually transforms to a ring rather than to an ellipsoid, and fragments at n c ~ 
7 x 10 18 cm~ 3 into two pieces with almost equal mass of ~ 1 x 10~ 2 M (Fig. [SH). The 
rapid drop of the normalized angular velocity at n c ~ 10 18 cm -3 (dashed line of Fig. [6]a) can 
be attributed to a rapid increase in the oblateness of the core (i.e., deformation to a flat 
disk shape), and consequent redistribution of the angular momentum of the core due to this 
non-axisymmetric density distribution. 

When the density reaches n ~ 10 20 cm -3 , a stellar core forms in each fragment. Each core 
is surrounded by a shock (Fig. [8]e and f). Although the stellar cores are initially ellipsoidal in 
shape (FigJSJe), they subsequently become spherical (Fig. [SJf). The mass of each stellar core 
reaches 3.1 x 10~ 3 M Q each at the end of the calculation. The two cores continue to part, 
and the final separation is R sep = 0.36 AU. Although the main bodies of the cores approach 
spherical in shape after the formation, their envelopes become more elongated, like tails. 
They eventually become detached (see the upper panel of Fig. [8]f). Having density peaks 
inside [at (x,y) ~ (±0.17, ±0.24)], these "tails" may hatch other stellar cores. 

The evolutionary sequences of number density and velocity before fragmentation are 
shown in Figures [9]a and b, respectively. Despite the deformation to a disk-like structure, 
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the evolution of the cor e is well described by the Larson- Periston self-similar solution until 



been found by iMatsumoto et al.l (119971 ) 



the fragmenta tion fc.f.. ISaigo et al.l 120001) . Similar behavior in the isothermal collapse has 



4-1.3. Bar-mode Fragmentation: Model 8 

Here, we see model 8, whose parameters are (a , A), A^) = (0.83, 10~ 2 , 0.3). The initial 
core has non-axisymmetric perturbation 30 times larger than in model 6, while the other 
parameters are the same. The evolutionary sequence of this model is shown in Figure [THl A 
non-axisymmetric structure, or a "bar" , already appears in an early phase. The bar becomes 
more elongated with collapse and eventually fragments. Evolution of the normalized angular 
velocity (u> c ), oblateness (e b), and axis ratio (e ar ) are plotted in Figures [TTb. b, and c (dashed 
lines) . 

Like in the previous examples, the central portion evolves keeping a spherical shape until 
disk formation (Fig. [Tol a; n c = 4.7 x 10 8 cm -3 ). Evolution of both the normalized angular 
velocity and the oblateness is also very similar to those of previous models (Fig. [6ja, b) . The 
oblateness saturates around e b — 7 at n c ~ 10 13 cm -3 and it subsequently oscillates (dashed 
lines of Fig. [TTb ) . When the oblateness reaches its maxima, a thin disk bounded by shocks 
forms, as in model 6 (middle and lower panels of Fig. [TUb ) . However, owing to its high initial 
value, the non-axisymmetric perturbation has grown significantly by the time of the disk 
formation, i.e., the saturation of uj c . The growth of the axis ratio is accelerated even more 
after the disk formation and a clear bar-like structure develops (Fig. [TUB) . The bar elongates 
further and eventually fragments at n c = 1.9 x 10 19 cm -3 . By this moment, the axis ratio 
reaches e ar = 31.4. This value is similar to the axis ratio at the fragmentation found by 
Tsuribe & Omukai (2006), whose calculation is limited to non-rotating cases. Inside the 
fragments, the stellar cores form after the equation of state becomes adiabatic (Fig. [TOl e. f). 

The separation between the fragments is R sep = 0.4 AU and the mass of a fragment 
is 3 x 10~ 2 M Q at the end of the calculation. The fragmentation epoch and separation are 
similar to those in model 6, which implies that the initial amplitude of the non-axisymmetric 
perturbation does not affect these quantities. In contract to the case of model 6, the frag- 
ments approach each other after fragmentation as a result of their small orbital angular 
momentum. Judging from the final state of our calculation, merger of the fragments appears 
inevitable although it has not yet occurred within our run. 
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4.2. Effect of rotation parameter fi Q 

In this section, we summarize effects of different initial rotation parameters (3$ on the 
evolution. As examples, in the second column of Figure [31 we compare six models 2, 6, 10, 
14, 18 and 22, which have different (3 (= 1(T 6 - 0.1) but with the same A^(= 10" 2 ). The 
normalized angular velocities, oblatenesses, and axis ratios are plotted in Figures EK W>, EEc 
for four of them (models 2, 6, 18 and 22). 

Except for the fastest rotation case of model 2, normalized angular velocities increase in 
the low density regime as u c oc n c , which is the spherical collapse relation. The normalized 
angular velocity saturates around uj c ~ 0.2 — 0.3 in all cases. Note that the saturation 
values of u c are similar in all models, although the epoch of saturation is earlier (i.e., at 
lower densities) for those with larger initial rotation. After the saturation, the normalized 
angular velocity oscillates around the saturation value u c ~ 0.2 — 0.3 until fragmen t ation . 



Sim ilar behaviors o f the an gular velocity have been found by iMatsumoto fc Hanawal (120031 ) 



and iMachida et al.l (j2006al ). who showed that the normalized angular velocity converges to 
a certain value in the self-similar collapse. In the case with the fastest rotation (model 2; 
/3o = 0.1), since the saturation level u c ~ 0.2 — 0.3 (/? ~ 0.04 — 0.09) has been reached from 
the onset, the normalized angular velocity remains almost constant until fragmentation. 

The shape of the cores remains almost spherical, namely, the oblateness is small e b — 1, 
until the saturation of the normalized angular velocity. It is only after the saturation of uj c 
that the oblateness begins to increase remarkably (Fig. [6]a, b). In the largest rotation case 
of model 2, the oblateness can increase from the onset since the saturation level of ui c has 
already been attained already in the initial state. Like u c , the oblateness also saturates at 
e oh ~ 10 and oscillates thereafter around this value until fragmentation. 

From this behavior of oblateness, we conclude that the saturation of u c to ~ 0.2 — 0.3 
is necessary for the formation of a thin (e ^ ~ 10) disk. A thin disk can continue collapsin g 



with conserved angular momentum for effective adiabatic index 7 < 1 (ISaigo et al.l 120001 ) . 
For 7 > 1, the central portion of the disk eventually supported by rotation and its collapse is 
halted, leading to fragmentation of the disk. In our case, where 7 > 1 in most of the density 
range, such a thin disk eventually fragments. Thus the saturation of u c is a prerequisite for 
fragmentation of the core. 

The evolution of the axis ratio is presented in Figure E]c. Also shown by a straight 
line is that for the isothermal gas in the linear theory e w oc nj (see Sec. 4.1.1). Before 
the formation of thin disks, the axis ratios decrease for n c < 10 8 cm -3 and then increase in 
higher densities by obeying the linear theory. After the flat disk with e ^ > 7 — 10 is formed, 
the growth rate of the axis ratio is more enhanced than the linear theory predicts. In fact, 



-12- 



the axis ratio continues to increase even when the polytropic index 7 is larger than 1.1 for 
n c > 10 16 cm" 3 . A similar effect is known in the isothermal collap se: a collapsing disk exhibits 



a larger growth rate of the ba r-mode perturbation than a sphere (JNakamura fc Hanawall 19971 ; 



Matsumoto fc Hanawalll999f ) 



Since fragmentation occurs only after the disk formation, cores with larger initial rota- 
tion tend to fragment at lower densities, and the separations between fragments are wider. 
For example, in model 2 (f3 = 0.1), with the fastest rotation, the fragmentation density is 
2 x 10 11 cm -3 and the separation between the fragments is R sep = 600 AU at the end of the 
calculation, while in model 18 (/3o = 10" 5 ), they are 6 x 10 20 cm -3 and 0.09AU, respectively. 
For slower rotation (j3 = 10~ 6 ; model 22), fragmentation does not take place. We also 
calculated a few cases with even smaller rotation energies and confirmed that the cores do 
not fragment. 



4.3. Effect of bar-mode perturbations 

Next, in order to study the effect of initial bar- mode perturbation A<p, we focus on four 
models (models 5-8) with the same (3 (= 10" 2 ) but with different A^ (= 0.001, 0.01, 0.1, 
and 0.2 for models 5, 6, 7, and 8, respectively). Their final states are presented in the second 
row of Figure [3 

With increasing A^, the mode of fragmentation changes from the ring to bar mode. 
In models with small A§ (< 0.01; models 5 and 6), a ring forms near the center owing to 
centrifugal force and subsequently fragments into several pieces. On the other hand, for 
larger A^ (> 0.1; models 7 and 8), the central region transforms from a disk to a bar-like 
structure before fragmentation. In model 7, the core takes a ring-like structure near the 
center whereas a bar-like one in a larger scale. The model with largest A^ (Model 8) clearly 
takes a bar-like structure. 

Whether fragmentation occurs or not appears to be solely determined by the rotation 
parameter (3q. Cores with (3q > 10~ 6 always experience fragmentation before the formation 
of stellar cores. In addition, fragmentation epoch and separation between fragments depend 
only weakly on bar- mode perturbation A^. The fragmentation density rif (i.e., epoch) ranges 
within an order of magnitude despite more than two orders of magnitude difference in A^\ 
m = 4 x 10 18 cm" 3 (A+ = 0.001; model 5), 7 x 10 18 cm" 3 (0.01; model 6), 7 x 10 18 cm' 3 (0.1; 
model 7), and 2 x 10 19 cm~ 3 (0.2; model 8). We can see a weak trend that for larger A^, 
the fragmentation densities are higher. The larger bar-mode perturbation A^ extracts more 
angular momentum from the central region, thereby delaying formation of a quasi-rotation- 
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supported disk and thus fragmentation. 

The evolution of normalized angular velocity u c , oblateness e b, and axis ratio e ar are 
shown in Figures [TTl a. b, c for models 5, 6, 7, and 8. With the same /3q, evolutionary tracks 
for both uo c and e Q b are very similar for all these models until n c ~ 10 15 cm~ 3 (Fig. [TTl a. 
b). For higher but still below fragmentation densities (n c ~ 10 15 — 10 18 cm -3 ), these tracks 
deviate from each other slightly owing to different shapes of the cores around the center 
(Fig. [3]): disks remain highly circular in low A^ models (5, 6 and 7), while an elongated bar 
forms in model 8 with high A^. As seen in Figure [TTJc, all the curves for £ ar have the same 
shape and are only shifted upward or downward in relation to one another. That is, their 
growth rates are common but initial values are different. Fragmentation occurs via the ring 
mode in models 5-7 while via the bar-mode in model 8. At the fragmentation, the axis ratios 
are less than unity for models 5-7 (e ar = 4.2 x 10 -2 , 0.47, and 0.46 for models 5, 6, and 7, 
respectively), whereas it reaches as high as 31.4 for model 8. 

In summary, the initial amplitude of bar-mode perturbation does affect the fragmenta- 
tion mode, but the fragmentation epoch and separation of fragments are mostly determined 
by the rotation parameter /3 Q . 



5. Evolution of Highly Unstable Cores: 
Effect of initial thermal energy ao 

In the previous section, we investigated the evolution of cores with thermal energy com- 
parable to gravitational energy (ao ~ 1). Although this setting is probably most realistic for 
initial cores according to numerical simulations of gravitational fragmentation (e.g., Bromm 
et al. 2002; Yoshida et al. 2006), to see dependence of core evolution on initial thermal 
energy parameter ao, we study cases with ao <C 1 in this section. 

Figure [T21 shows the core structure at the end of the calculation for models 25-28. These 
models have the same initial rotation parameters (j3q = 10~ 3 ) and bar-mode perturbations 
{Atj, = 0.01), but different thermal energies (ao)- Figure [131 shows the evolution of normalized 
angular velocity u> c , oblateness e b, and axis ratio £ ar . The difference in «o hardly affects the 
evolution of u c and e b, although, for models with smaller a , both the u c and e b saturate 
at slightly higher densities and at slightly higher values (Fig. [T3"l a. b). On the other hand, 
the axis ratios e ar evolve quite differently among these models, in particular, in low densities 
n c < 10 10 cm -3 (Fig. [T3lc): e ar increases more for the lower ao cores. For higher densities, 
where the evolution of the cores already converge to the self-similar solution, the growth rate 
of e ar becomes the same for all the models, but this convergence occurs later for the lower 
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«o cores. Therefore, cores with smaller olq tend to have larger e ar , and then to fragment in 
the bar mode. For example, in models 25 and 26 with small a^ (=0.2, 0.4), the bar-mode 
perturbation grows sufficiently and the cores fragment in the bar mode (Fig. [T2l a. b), while 
in models 27 and 28 (ao = 0.6 and 0.8), fragmentation occurs in the ring mode (Fig. [T2l c. 
d). 



6. Fragmentation Condition 

As seen in Figure EJ the condition for fragmentation appears to be 

A, > /3crit ~ 10- 6 - 10- 5 , (8) 

i.e., fragmentation occurs if the initial core has rotation energy 10~ 6 -10 -5 times the grav- 
itational energy. Or, in terms of angular velocity of the initial core, this can be rewritten 



fio > 3.67 x 10- 17 (— ^S) 2/ V\ (9) 

\10 d cm - ' 5 / 



as 

_^o ^ 3 

.10 3 cm -3 - 

where Qq, and no are the angular velocity and the central number density of the initial 

core, respectively, and we set j3 crit = 10~ 6 . The dependence of equation (jUJ) on the density 

comes from the fact that the angular velocity increases Q oc n c in the case of a spherically 

collapsing core. 

Next, we present a physical explanation for the fragmentation condition (jHJ). As in 
Sec. 4. 2 and 4.3, the formation of a thin disk with oblateness e b — 10 is necessary for 
fragmentation. This must be attained before the equation of state becomes adiabatic at 
n c ~ 10 20 cm -3 since the core becomes spherical and fragmentation is prohibited thereafter. 
The oblateness begins to grow remarkably only after the normalized angular velocity becomes 
close to the saturation value ~ 0.2 — 0.3, and the level of e b — 10 is reached after further 
contraction of about four or five orders of magnitude in density. Therefore, as a condition 
for fragmentation, the normalized angular velocity must reach the saturation value (u c ~ 
0.2 — 0.3, or (3 C ~ 0.1) before ~ 10 15 — 10 16 cm -3 . By using the relation j3 c = y/uJ^ oc nj for 
the spherical collapse, this can be translated to that on initial (i.e., at ~ 10 3 cm~ 3 ) rotation 
parameter as /3q > (0.5 — 1) X 10~ 5 . This condition coincides well with our empirical relation 
©. 

In deriving the above condition, we neglect angular momentum transfer before the 
formation of a thin disk. As seen in Sec. 4.3, with the higher the bar-mode perturbation 
Afa the more angular momentum is transported. Thus, cores with higher initial A^ would 
need higher (3q for fragmentation to compensate the angular momentum transfer. This effect 
appears to be modest and is not observed clearly in our calculation. 
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7. Summary and Discussion 

We have studied the gravitational collapse of zero-metallicity star-forming cores for 
different initial conditions, which are specified by the following three parameters, i.e., thermal 
energy, rotation energy, and amplitude of non-axisymmetric perturbations. Since the thermal 
energy of star-forming cores is expected to be similar to the gravitational energy, particular 
attention is paid to such cores. Initial conditions of the cores where fragmentation occurs are 
surveyed. We found that the cores with initial rotation energy higher than 10 -6 — 10~ 5 times 
the gravitational energy fragment into a few pieces. Dependence on the other parameters 
is weak and is not observed clearly in our calculation. The physical explanation for this 
fragmentation condition is also presented. 

How much angular momentum d o initial cores have in reality? According to cosmological 



simulations by lYoshida et al.l (120061 ). first-star forming cores have the rotation parameter 
j3 ~ 0.1. In discussing the fragmentation condition, the rotation parameter near the center 
should be used rather than that averaged over the core. Matter in the outer part tends to 
have larger rotation energy than gravitational energy, the value of the rotation parameter j3 
becomes larger if outer material is included in calculations. Even though the central value 
of the rotation parameter might be significantly smaller than 0.1 owing to this effect, it is 
probably higher than our threshold value for fragmentation (10~ 5 — 10~ 6 ) with a wide margin. 
Therefore, most of the cores are expected to fragment and proto-binaries are formed even in 
primordial-gas clouds. 

In Figure[EH we plot the separations between fragments R sep at the end of the calculation 
against the fragmentation epoch n^ for all the models where fragmentation is observed. The 
large crosses show ranges of fragmentation epochs (horizontal) and fragmentation scales 
(vertical) for models with (3 > 0.02 (dotted line) (j3o < 0.02, solid line, respectively). 
The Jeans length (left axis) and gas temperature (right axis) are also shown. We can see 
that the separation between fragments is approximately equal to the Jeans length at the 
fragmentation. Another striking feature is that fragmentation densities are clustered in two 
epochs of n « 10 12 cm -3 and 10 21 cm -3 . Cores with rapid rotation (/? > 0.02; symbols '+') 
fragment in the earlier epoch (n c ~ 10 10 — 10 14 cm" 3 ) with larger separation ~ 10 — 1000 AU, 
while those with slower rotation (/3q < 0.02; 'x') in the later epoch (n c ~ 10 18 — 10 21 cm~ 3 ) 
with smaller separation ~ 0.01 — 1AU. These two fragmentation epochs have their origin 
in thermal evolution of primordial gas. The lower-density fragmentation is a result of a 
temperature drop owing to the three-body H2 formation at n c ~ 10 10 cm -3 , while the higher- 
density one results from another temperature drop at n c ~ 10 17 cm -3 , owing to the H2 
collision-induced cooling. Recall that a quasi-rotation supported disk can continue to collapse 
for 7 < 1 while conserving angular momentum. Thus even if collapse of the disk is about 
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to halt, the disk can collapse further owing to the temperature drops. Such disks stop their 
collapse and fragment when the temperature starts increasing again. It is known that these 
two cooling mechanisms do not cause fragmentation of the core by simple thermal instability 
(Omukai & Yoshii 2003; Ripamonti & Abel 2004). In our case, where a disk is present, these 
cooling mechanisms result in the fragmentation of the disk. 

There is a caveat in our results on fragmentation epochs. We adopted the equation of 
state derived from one-zone calculation. For example, the temperature variation in one- zone 
models tends to be exaggerated: the two temperature dips are deeper than in hydrody- 
namical calculations (e.g., Omukai & Nishi 1998). In addition, in hydrodynamical models, 
temp erature in the envelope tends to be higher than t hat in the center for the same den- 
sity flWhitehouse &: Batdl2006l ; IStamatellos et al.l 120071 ). Despite these possible differences, 
we expect that the fragmentation condition itself (eq. ED remains valid. In deriving this 
condition (Sec. 6), we only assumed that once a quasi- rotation supported disk forms, it 
eventually fragments. This is correct if the EOS has a sufficiently large density range with 
7 > 1, which is indeed the case in the primordial gas. In fact, in the context of present-day 
star formation, the same fra gmentation criterion has been reported ( jMatsumoto &: Hanawa 
20031 ; iMachida et al.ll2005bl ). On the other hand, the fragmentation epoch and scale can be 
subject to change if EOS is altered. 

Although the majority of cores are predicted to fragment before the stellar core for- 
mation, whether they finally evolve into bona fide binaries/multiples remains obscure. The 
separation between fragments changes significantly during the mass accretion phase as a 
result of angular momentum influx to the fragments. An accreted gas in a later stage tends 
to have a larger specific angular momentum, while spiral arms in a density pattern, which 
transfer angular momentum outward, might be invoked in such a phase. To answer whether 
fragments evolve into binaries/multiples or merge into a single star, we need to extend cal- 
culations over the entire accretion phase. 

When more than two fragments appear (e.g., models 1 and 2 in Fig. [3]), some m i ght b e 
ejected from the parental core as a result of three-body interactions (JBate et al.ll2002l . 120031 ). 
Ejected objects stop gaining mass by accretion thereafter. By this mechanism, a spectrum 
in mass of metal-free stars might be formed. For example, if an ejection event takes place 
in the very early accretion phase, the ejected object can eventually evolve into a metal-free 
brown dwarf. Or if a stellar core with its envelope of ~ 1 M & is ejected, a low-mass metal-free 
star is the outcome. Considering the origin of the two hyper m etal-poor ([Fe/H]< —5) stars 
(Christlieb et al. 2001; Frebel et al. 2005). ISuda et al.l (120041 ) showed that their abundance 
patterns can be explained by nucleosynthesis and mass transfer in first-generation low-mass 
binary stars. In our scenario, such a system can be formed naturally. In some models (e.g., 
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see Fig. [3]a and b), multiple pairs of binaries appear. If one such pair is ejected from the 
parental core, it would result in a low-mass metal-free binary system. 

With neither ejection nor merging to a single star, our binary systems of fragments are 
expected to become very massive star binaries. Such systems are considered to be important 
sources of gravitational waves (Belczy nski et al. 2004). Also they can be gamm a-ray burst 
progenitors if formed in close binaries (iBromm fc Loeb 1 120061 ; lYoon et al.l 120061 ) . 



In this study, we did not include magnetic effects. According to rece nt studies, however 



a low level of mag netic field could be present even in the early universe (JLanger et al.l 12003 
Ichiki et al.ll2006l ). If an initial primordial-gas core h as B > 10~ 13 ' 5 G of the magnetic field 



a protostellar jet is launched (IMachida et al.l 12006a ). and this jet significantly affects the 
mass accretion process. Even with a lower level, a magnetic field can suppress fragmentation 



of the core, by transferr ing angular momentum via magnetic braking (e.g., IMachida et al. 
20071 ; iPrice fc Batdl2007l ). So far those effects on primordial star formation have been largely 
overlooked. Further studies along this line are clearly needed in future. 
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Table 1: Model parameters and calculation results 
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furthermost fragments, d fragmentation mode (R: ring, B: bar). 
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Fig. 1.— ■ Gas pressure normalized by number density derived from a one-zone calculation 
(solid) and our fit used in the numerical simulation (dashed line). 
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Fig. 2. — The density (thick solid line) and enclosed mass (thick dashed line) distribution 
for the critical Bonnor-Ebert sphere as a function of the radius. Those increased by a factor 
/=1.4 are also plotted by the thin lines. 




Fig. 3. — The collapse outcome of primordial-gas cores with different initial rotation param- 
eter /3 and amplitude of non-axisymmetric (bar-mode) perturbation A^. The final states in 
our calculation for models 1-24 in Table 1, where the cores are close to hydrostatic initially, 
i.e., «o ~ 1 are presented. Shown by false colors is the density distribution around the 
center on the plane perpendicular to the rotation axis. The model numbers are denoted on 
the upper-left corners of the panels. The name of the panel (a-x) , the time elapsed since the 
beginning of the calculation t, the maximum number density attained n max are given at the 
top, and the box scale is shown at the bottom of the panel. 
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Fig. 4. — A schematic for the evolution of rotating cores. If the core reaches the stellar 
density p* before the disk formation, it does not fragment. Otherwise, it fragments in the 
bar mode or ring mode depending on the amplitude of the bar-mode perturbation at the 
time of disk formation. 
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Fig. 5.— The density {color scale and contour) and velocity (arrows) distributions at six 
different epochs in model 22 [(a , Po, A<f,) = (0.83, 10~ 6 , 10 -2 )]. The upper, middle, and 
lower panels show those on x — y, x — z, and y — z planes, respectively. The elapsed time, 
maximum density (above) and arrow scale (below) are denoted above or below each panel. 
The level of the grid is indicated in the upper left corner of each panel. Solid squares in each 
panel are outer boundaries of subgrids. The structure around the center is zoomed up from 
Panel a toward /. The grid size of Panel a (the coarsest grid shown is I = 2) is 6.7 x 10 6 AU, 
while that of Panel / ' (I = 27) is 0.2 AU, 2 25 times finer than in Panel a. The black contour in 
Panel /denotes the constant density surface of n — 10 20 cm -3 , approximately corresponding 
to the surface of the stellar core. 
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Fig. 6. — Evolution of (a) the central normalized angular velocity u c , (6) oblateness e b 
and (c) axis ratio e ar of the cores as a function of the central number density n c . Four 
cases (models 2, 6, 18 and 22 in the second column of Figure [3]) with the same initial 
non-axisymmetric perturbation A^ = 0.01 but with different rotation parameters (/3q = 
10~ 6 — 10 _1 ) are shown. The dash in Panel (a) shows the growth rate for the spherical 
collapse, while that in Panel (c) indicates the growth rate in the isothermal case by the 
linear analysis. 
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Fig. 7.— The radially averaged (a) number density and (b) velocity profiles at different 
epochs in model 22. The time for each plot is indicated in panel a. The dash in panel a 
shows the relation n oc r -2 ' 2 . 
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Fig. 8.— The same as Fig. but for model 6 with (a , A), A^) = (0.83, 10~ 2 , 0.01). 
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Fig. 9.— The same as Figure [7] but for model 6. 
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Fig. 10.— The same as Fig. El but for model 8 with (a , Po, A$) = (0.83, 10~ 2 , 0.3). 
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Fig. 11. — The same as Figure[6]but for cases with same initial rotation parameter /3q = 10~ 2 
with different amplitude of initial non-axisymmetric perturbation A$ = 0.01 (models 5, 6, 7 
and 8, in the second row of Figure [3]). 
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Fig. 12.— The collapse outcome of highly gravitationally unstable cores (models 25-28). 
The values of «o are indicated in the bottom-left corner of the panels. The legend is the 
same as in Figure El 
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Fig. 13. — The same as Figure [6]but for highly gravitationally unstable cores (models 25-25 
The values of «o are indicated in panel (b). 
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Fig. 14.— The separations between fragments against the density at fragmentation. We 
adopt the separation between furthermost fragments if more than two fragments appeared. 
The evolution of temperature (right) and Jeans length (left) is also illustrated by the dotted 
and solid lines, respectively. The results of all the models where fragmentation is observed 
are plotted. Those with the rotation parameter /3 > 0.02 (< 0.02) are shown by symbols 
'+', ('x', respectively). The large dotted (solid) cross represents the range of the fragment 
separation and fragmentation epoch for cases with (3q > 0.02 (< 0.02, respectively). 



